

## save stata data file (Rdata file not available at time) to working directory. Here I preserved the name UnitedStates LAPOP AmericasBarometer 2019 v1.0_W.dta

#read file
install.packages('readstata13')
library('readstata13')
abus19<-read.dta13("UnitedStates LAPOP AmericasBarometer 2019 v1.0_W.dta", convert.factors=F)


install.packages('survey')
library(survey)

install.packages('jtools')
library(jtools)

install.packages('weights')
library(weights)



#demographics

abus19$sex<-abus19$q1-1

abus19$ageyears<-abus19$q2


abus19$hhincome<-abus19$q10new/16

abus19$collegedegree[abus19$ed_usa %in% c(5, 6)]<-1
abus19$collegedegree[!(abus19$ed_usa %in% c(5, 6))]<-0


abus19$blackvsall[abus19$etid==4]<-1
abus19$blackvsall[abus19$etid!=4|is.na(abus19$etid)]<-0
abus19$latinovsall[abus19$etid==4003]<-1
abus19$latinovsall[abus19$etid!=4003|is.na(abus19$etid)]<-0
abus19$asianvsall[abus19$etid==4004]<-1
abus19$asianvsall[abus19$etid!=4004|is.na(abus19$etid)]<-0
abus19$nativeamvsall[abus19$etid==4005]<-1
abus19$nativeamvsall[abus19$etid!=4005|is.na(abus19$etid)]<-0
abus19$mixedvsall[abus19$etid==4006]<-1
abus19$mixedvsall[abus19$etid!=4006|is.na(abus19$etid)]<-0
abus19$middleeastvsall[abus19$etid==4008]<-1
abus19$middleeastvsall[abus19$etid!=4008|is.na(abus19$etid)]<-0
abus19$othervsall[abus19$etid==7|is.na(abus19$etid)]<-1
abus19$othervsall[abus19$etid!=7]<-0

#openness to authoritarian governance items

#half-sample only
abus19$closecong[abus19$jc15a==1]<-1
abus19$closecong[abus19$jc15a==2]<-0

#half-sample only
abus19$closesupsup[abus19$jc16a==1]<-1
abus19$closesupsup[abus19$jc16a==2]<-0

abus19$closeinstit<-rowMeans(data.frame(abus19$closecong, abus19$closesup), na.rm=T)




abus19$churchillrev<-1-((abus19$ing4-1)/6)


abus19$disallowvote<-1-(abus19$d1-1)/9
abus19$disallowdemonstration<-1-(abus19$d2-1)/9
abus19$disallowrunforoffice<-1-(abus19$d3-1)/9
abus19$disallowspeeches<-1-(abus19$d4-1)/9


#economic left items

abus19$incomeineq<-(abus19$ros4-1)/6

abus19$helppoor<-(abus19$redist1-1)/6

abus19$richtax<-(abus19$redist2a-1)/6

abus19$unemp<-1-((abus19$redist3-1)/6)



#cultural conservatism items

abus19$religattend<-1-((abus19$q5a-1)/4)

abus19$religimport<-1-((abus19$q5b-1)/3)

abus19$menbetterleaders<-1-((abus19$vb50-1)/3)





#for control variables


abus19$trustcongress<-(abus19$b13-1)/6
abus19$trustparties<-(abus19$b21-1)/6
abus19$trustelections<-(abus19$b47a-1)/6



abus19$satisdemoc<-1-((abus19$pn4-1)/3)


abus19$polint<-1-((abus19$pol1-1)/3)


abus19$follownews<-1-((abus19$gi0-1)/4)



#descriptives


wtd.mean(abus19$churchillrev, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$churchillrev, weights=abus19$wt, na.rm=T))

abus19$intol<-with(abus19, rowMeans(data.frame(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches)), na.rm=T)

wtd.mean(abus19$intol, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$intol, weights=abus19$wt, na.rm=T))

abus19intolitems<-subset(abus19, select=c(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches))

psych::alpha(abus19intolitems)


wtd.mean(abus19$closeinstit, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$closeinstit, weights=abus19$wt, na.rm=T))



abus19$authgov<-with(abus19, rowMeans(data.frame(intol, closeinstit, churchillrev), na.rm=T))

wtd.mean(abus19$authgov, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$authgov, weights=abus19$wt, na.rm=T))

abus19authgovitems<-subset(abus19, select=c(intol, closeinstit, churchillrev))

psych::alpha(abus19authgovitems)



abus19$cultright<-with(abus19, rowMeans(data.frame(religattend, religimport, menbetterleaders), na.rm=T))



abus19$cultrightnoattend<-with(abus19, rowMeans(data.frame(religimport, menbetterleaders), na.rm=T))

wtd.mean(abus19$cultright, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$cultright, weights=abus19$wt, na.rm=T))

abus19cultrightitems<-subset(abus19, select=c(religattend, religimport, menbetterleaders))

psych::alpha(abus19cultrightitems)



abus19$econleft<-with(abus19, rowMeans(data.frame(incomeineq, helppoor, richtax, unemp), na.rm=T))

wtd.mean(abus19$econleft, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$econleft, weights=abus19$wt, na.rm=T))

abus19econleftitems<-subset(abus19, select=c(incomeineq, helppoor, richtax, unemp))

psych::alpha(abus19econleftitems)



abus19$confinstit<-with(abus19, rowMeans(data.frame(trustelections, trustparties, trustcongress)),na.rm=T)

wtd.mean(abus19$confinstit, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$confinstit, weights=abus19$wt, na.rm=T))

abus19confinstititems<-subset(abus19, select=c(trustelections, trustparties, trustcongress))

psych::alpha(abus19confinstititems)



abus19$poleng<-with(abus19, rowMeans(data.frame(polint, follownews)),na.rm=T)

wtd.mean(abus19$poleng, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$poleng, weights=abus19$wt, na.rm=T))

abus19polengitems<-subset(abus19, select=c(polint, follownews))

psych::alpha(abus19polengitems)



wtd.mean(abus19$satisdemoc, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$satisdemoc, weights=abus19$wt, na.rm=T))



wtd.mean(abus19$collegedegree, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$collegedegree, weights=abus19$wt, na.rm=T))




wtd.mean(abus19$ageyears, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$ageyears, weights=abus19$wt, na.rm=T))




wtd.mean(abus19$sex, weights=abus19$wt, na.rm=T)

sqrt(wtd.var(abus19$sex, weights=abus19$wt, na.rm=T))


median(abus19$q10new, na.rm=T)





#main analyses




library(mice)
library(miceadds)
library(survey)

xtract.svymi <- function(z){
  #require(norm)
  m2 <- length(z)
  ns <- nobs(z[[1]]) #first dataset
  betas <- lapply(z, FUN = function(x){coef(x)} )
  sess <- lapply(z, FUN = function(x){vcov(x)}) #the vcov
  #ses <- lapply(sess, FUN = function(x){sqrt(diag(x))}) #the se, not the complete vcov
  
  dm <- lapply(z, FUN = function(x){model.matrix(x)})
  r.1 <- sapply(1:m2, function(row) cor(dm[[row]] %*% betas[[row]], z[[row]]$y))
  fisher.r2z <- function(r) {0.5 * (log(1+r) - log(1-r))}
  zz <- mean(fisher.r2z(r.1))
  r2 <- psych::fisherz2r(zz)^2 #in psych
  R2 <- mean(r2)
  
  test1 <- summary(pool_mi(betas, sess))
  
  test1$term <- row.names(test1)
  test2 <- test1[,c(8,1,2,4)] #4 is the pvalue
  
  row.names(test2) <- c()
  names(test2) <- c('term','estimate','std.error', 'p.value')
  
  tr <- createTexreg(
    coef.names = test2$term, 
    coef = test2$estimate, 
    se = test2$std.error, 
    pvalues = test2$p.value,
    gof.names = c("R2", "Nobs"), 
    gof = c(R2, ns),
    gof.decimal = c(T,F)
  )
  
}


abus19<-subset(abus19, select=c(sex, ageyears, collegedegree, hhincome, blackvsall, latinovsall, asianvsall, 
                            nativeamvsall, mixedvsall, middleeastvsall, othervsall, closeinstit, churchillrev, disallowvote, 
                            disallowdemonstration, disallowrunforoffice, 
                            disallowspeeches, incomeineq, helppoor, richtax, unemp, religattend, religimport, menbetterleaders, 
                            trustelections, trustparties, trustcongress, polint, follownews, satisdemoc, estratopri, wt, upm))


abus19i<-mice(abus19, m=20)

longabus19 <- complete(abus19i, action='long', include=TRUE)

longabus19$intol<-with(longabus19, rowMeans(data.frame(disallowvote, disallowdemonstration, disallowrunforoffice, disallowspeeches)))

longabus19$authgov<-  with(longabus19, rowMeans(data.frame(closeinstit, churchillrev, intol)))


longabus19$cultright<-  with(longabus19, rowMeans(data.frame(religattend, religimport, menbetterleaders)))


longabus19$econleft<-  with(longabus19, rowMeans(data.frame(incomeineq, helppoor, richtax, unemp)))

longabus19$age1s<-with(longabus19, (ageyears-min(ageyears, na.rm=T))/(max(ageyears, na.rm=T)-min(ageyears, na.rm=T)) )

longabus19$confinstit<-with(longabus19, rowMeans(data.frame(trustelections, trustparties, trustcongress)))

longabus19$poleng<-with(longabus19, rowMeans(data.frame(polint, follownews)))


longabus19$cultrightplussd<-with(longabus19, cultright-(mean(cultright, na.rm=T)+sd(cultright, na.rm=T)))

longabus19$cultrightminussd<-with(longabus19, cultright-(mean(cultright, na.rm=T)-sd(cultright, na.rm=T)))

longabus19$econleftplussd<-with(longabus19, econleft-(mean(econleft, na.rm=T)+sd(econleft, na.rm=T)))

longabus19$econleftminussd<-with(longabus19, econleft-(mean(econleft, na.rm=T)-sd(econleft, na.rm=T)))


longabus19$cultrightscl <- with(longabus19, (cultright-mean(cultright, na.rm=T))/(2*sd(cultright, na.rm=T)))

longabus19$econleftscl <- with(longabus19, (econleft-mean(econleft, na.rm=T))/(2*sd(econleft, na.rm=T)))

longabus19$polengscl <-  with(longabus19, (poleng-mean(poleng, na.rm=T))/(2*sd(poleng, na.rm=T)))


longabus19$age1sscl <- with(longabus19, (age1s-mean(age1s, na.rm=T))/(2*sd(age1s, na.rm=T)))


longabus19$satisdemocscl <- with(longabus19, (satisdemoc-mean(satisdemoc, na.rm=T))/(2*sd(satisdemoc, na.rm=T)))


longabus19$confinstitscl <- with(longabus19, (confinstit-mean(confinstit, na.rm=T))/(2*sd(confinstit, na.rm=T)))


longabus19$hhincomescl <- with(longabus19, (hhincome-mean(hhincome, na.rm=T))/(2*sd(hhincome, na.rm=T)))


longabus19$econextrem<-abs(with(longabus19, (econleft-mean(econleft, na.rm=T))))

longabus19$econextremscl<-with(longabus19, (econextrem-mean(econextrem, na.rm=T))/(2*sd(econextrem, na.rm=T)))



abus19i2 <- as.mids(longabus19)


abus19_imp_list <- imputationList(lapply(1:abus19i2$m, function(n) mice::complete(abus19i2, action = n)))

abus19design<-svydesign(~upm,strata=~estratopri, nest=TRUE, fpc=NULL, data=abus19_imp_list, weights=~wt)


abus19reg1<- with(abus19design,svyglm(authgov~cultrightscl+econleftscl))


abus19reg2<- with(abus19design,svyglm(authgov~cultrightscl+econleftscl+cultrightscl:econleftscl))


abus19reg3<- with(abus19design,svyglm(authgov~cultrightscl+econleftscl+sex+age1sscl+collegedegree+hhincomescl+blackvsall+latinovsall+asianvsall+nativeamvsall+mixedvsall+middleeastvsall+othervsall+confinstitscl+polengscl+satisdemocscl))


abus19reg4<- with(abus19design,svyglm(authgov~cultrightscl+econleftscl+sex+age1sscl+collegedegree+hhincomescl+blackvsall+latinovsall+asianvsall+nativeamvsall+mixedvsall+middleeastvsall+othervsall+confinstitscl+polengscl+satisdemocscl+cultrightscl:econleftscl))




abus19regright<- with(abus19design,svyglm(authgov~cultrightplussd+econleftminussd+cultrightplussd:econleftminussd))


abus19regleft<- with(abus19design,svyglm(authgov~cultrightminussd+econleftplussd+cultrightminussd:econleftplussd))


abus19regprot<- with(abus19design,svyglm(authgov~cultrightplussd+econleftplussd+cultrightplussd:econleftplussd))


abus19regfree<- with(abus19design,svyglm(authgov~cultrightminussd+econleftminussd+cultrightminussd:econleftminussd))

abus19out1 <- xtract.svymi(abus19reg1)
abus19out2 <- xtract.svymi(abus19reg2)
abus19out3 <- xtract.svymi(abus19reg3)
abus19out4 <- xtract.svymi(abus19reg4)


htmlreg(list(abus19out1,abus19out2,abus19out3,abus19out4), digits = 3,file = "lapopUS19.html",custom.coef.names = c("Constant","Cultural Conservatism","Left Economic Attitudes","Cultural Conservatism x Left Economic Attitudes","Female","Age","College Degree","Household Income","Black","Latino","Asian","Native American","Mixed Race","Middle Eastern", "Other Race","Confidence in Institutions","Political Engagement","Satisfaction with Democracy"),caption.above = T,caption = "Table D-17: Predictors of Openness to Authoritarian Governance in LAPOP U.S.-2019 Dataset",custom.note = ("%stars. <br> Continuous predictors are mean-centered and scaled by two standard deviations. Binary predictors are coded 0-1.<br>The outcome variable is coded to range from 0.0 to 1.0 <br> Missing data were replaced using multiple imputation with results pooled across 20 copies of the dataset. <br> Comparison Racial-Ethnic Group is White"),custom.gof.names = c("R<sup>2</sup>","Observations"))

abus19cultest <- summary(pool(abus19reg1),conf.int = T)[2,2]
abus19cultlower <- summary(pool(abus19reg1),conf.int = T)[2,7]
abus19cultupper <- summary(pool(abus19reg1),conf.int = T)[2,8]

abus19econest <- summary(pool(abus19reg1),conf.int = T)[3,2]
abus19econlower <- summary(pool(abus19reg1),conf.int = T)[3,7]
abus19econupper <- summary(pool(abus19reg1),conf.int = T)[3,8]

abus19intest <- summary(pool(abus19reg2),conf.int = T)[4,2]
abus19intlower <- summary(pool(abus19reg2),conf.int = T)[4,7]
abus19intupper <- summary(pool(abus19reg2),conf.int = T)[4,8]


abus19rightest <- summary(pool(abus19regright),conf.int = T)[1,2]
abus19rightlower <- summary(pool(abus19regright),conf.int = T)[1,7]
abus19rightupper <- summary(pool(abus19regright),conf.int = T)[1,8]


abus19leftest <- summary(pool(abus19regleft),conf.int = T)[1,2]
abus19leftlower <- summary(pool(abus19regleft),conf.int = T)[1,7]
abus19leftupper <- summary(pool(abus19regleft),conf.int = T)[1,8]


abus19protest <- summary(pool(abus19regprot),conf.int = T)[1,2]
abus19protlower <- summary(pool(abus19regprot),conf.int = T)[1,7]
abus19protupper <- summary(pool(abus19regprot),conf.int = T)[1,8]


abus19freeest <- summary(pool(abus19regfree),conf.int = T)[1,2]
abus19freelower <- summary(pool(abus19regfree),conf.int = T)[1,7]
abus19freeupper <- summary(pool(abus19regfree),conf.int = T)[1,8]


abus19cultb3<-extract.cis(abus19reg3)[2,1]
abus19polengb3<-extract.cis(abus19reg3)[16,1]
abus19ageb3<-extract.cis(abus19reg3)[5,1]
abus19collegeb3<-extract.cis(abus19reg3)[6,1]
abus19confinstitb3<-extract.cis(abus19reg3)[15,1]



#additional analyses without imputation, scales with all avail items and listwise deletion in analysis.  Reported in Part E of SOM.


abus19$cultrightscl<-with(abus19, (cultright-mean(abus19$cultright, na.rm=T))/(2*sd(cultright, na.rm=T)))

abus19$econleftscl<-with(abus19, (econleft-mean(abus19$econleft, na.rm=T))/(2*sd(econleft, na.rm=T)))


abus19$cultrightplussd<-with(abus19, cultright-(mean(cultright, na.rm=T)+sd(cultright, na.rm=T)))

abus19$cultrightminussd<-with(abus19, cultright-(mean(cultright, na.rm=T)-sd(cultright, na.rm=T)))

abus19$econleftplussd<-with(abus19, econleft-(mean(econleft, na.rm=T)+sd(econleft, na.rm=T)))

abus19$econleftminussd<-with(abus19, econleft-(mean(econleft, na.rm=T)-sd(econleft, na.rm=T)))


abus19design<-svydesign(~upm, strata=abus19$estratopri, nest=TRUE, fpc=NULL, data=abus19, weights=abus19$wt)


reg1abus19l<- svyglm(authgov~cultrightscl+econleftscl, design=abus19design)


reg2abus19l<- svyglm(authgov~cultrightscl+econleftscl+cultrightscl:econleftscl, design=abus19design)

abus19lregright<- svyglm(authgov~cultrightplussd+econleftminussd+cultrightplussd:econleftminussd, design=abus19design)


abus19lregleft<- svyglm(authgov~cultrightminussd+econleftplussd+cultrightminussd:econleftplussd, design=abus19design)


abus19lregprot<- svyglm(authgov~cultrightplussd+econleftplussd+cultrightplussd:econleftplussd, design=abus19design)


abus19lregfree<- svyglm(authgov~cultrightminussd+econleftminussd+cultrightminussd:econleftminussd, design=abus19design)



abus19lcultest<-summary(reg1abus19l)$coef[2,1]
abus19lcultlower<-confint(reg1abus19l)[2,1]
abus19lcultupper<-confint(reg1abus19l)[2,2]

abus19leconest<-summary(reg1abus19l)$coef[3,1]
abus19leconlower<-confint(reg1abus19l)[3,1]
abus19leconupper<-confint(reg1abus19l)[3,2]

abus19lintest<-summary(reg2abus19l)$coef[4,1]
abus19lintlower<-confint(reg2abus19l)[4,1]
abus19lintupper<-confint(reg2abus19l)[4,2]


abus19lrightest<-summary(abus19lregright)$coef[1,1]
abus19lrightlower<-confint(abus19lregright)[1,1]
abus19lrightupper<-confint(abus19lregright)[1,2]

abus19lleftest<-summary(abus19lregleft)$coef[1,1]
abus19lleftlower<-confint(abus19lregleft)[1,1]
abus19lleftupper<-confint(abus19lregleft)[1,2]

abus19lprotest<-summary(abus19lregprot)$coef[1,1]
abus19lprotlower<-confint(abus19lregprot)[1,1]
abus19lprotupper<-confint(abus19lregprot)[1,2]

abus19lfreeest<-summary(abus19lregfree)$coef[1,1]
abus19lfreelower<-confint(abus19lregfree)[1,1]
abus19lfreeupper<-confint(abus19lregfree)[1,2]



